#!/usr/bin/perl -w
use strict;
use warnings;
use lib "/net/cpp-group/Leo/bin";
use parse_bl2seq;
use parse_fasta;
use Getopt::Long;

#originally
#parse_codeml_positions.pl		Scott Beatson 		20.05.03
# completely rewritten by Leo Goodstadt
#
#A program to parse the output from multiple codeml mlc files to generate
#a fasta formatted file that can be used as a footer in an alignment program

my $usage = <<'USAGE';

USAGE:

parse_codeml_results.pl site_specific_results_files
				--parameter_output output.parameters
				--positive_residues output.fa
				[--verbose]
				[--help]

    Parses the output from multiple codeml result files to generate
    a fasta formatted file indicating +ve selected positions.
	Codeml parameters are saved in the <parameter_output>

USAGE

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#   Get options

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
# mandatory


# optional parameters
my $help = undef;
my $verbose = undef;
my $positive_residues;
my $parameter_output;
GetOptions(
			'help'			   		=> \$help,
			'verbose'		   		=> \$verbose,
			'positive_residues=s'	=> \$positive_residues,
			'parameter_output=s'	=> \$parameter_output,
			);

die $usage if ($help);
die $usage unless ($parameter_output  && $positive_residues);
die $usage unless (defined($ARGV[0]));

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#	Get gaps in original sequences (removed by codeml)
#
=pod
my @common_seq;
sub get_seq_callback($$$$$$)
{
    my ($acc,  $seq, $save_sequences) = @_[0, 1, 5];
	
	my $sequence = join "", @$seq;

	if (!@common_seq)
	{
		@common_seq = ('.') x length($sequence);
	}
	die "Error:\n\tSequences difference length ".
			"(". @common_seq. ",".length($sequence).")" unless
		(@common_seq == length($sequence));

	my @aa = split //, $sequence;
	for (my $i = 0; $i < @aa; ++$i)
	{
		$common_seq[$i] = '-' if ($aa[$i] eq '-');
	}
}

open (PROT_ALIGNED, $prot_aligned) or
	die "Error\n\tCould not open the protein alignment FASTA file [prot_aligned]\n$!\n";
my $cnt_seq = parse_fasta::parse_sequences(  *PROT_ALIGNED,
												\&get_seq_callback,
												undef);
print STDERR "\n\t$cnt_seq\tsequences parsed.\n\n" if $verbose;
my $common_seq = join ("", @common_seq);
my @alignment = get_alignment_marks($common_seq);
=cut

#print STDERR gap_positions_array_to_str(@alignment), "\n";
#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888






#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#	test_model
#
#		make sure prior parameters of model also indicate +ve selection
#		i.e. chi significance test and that there is a +ve selected category
#
sub test_model($\%$$)
{
	my ($output_file, $data, $model1, $model2) = @_;
	die unless exists $data->{$model2};
	die unless exists $data->{$model1};
	die unless exists $data->{$model2}{'cnt_free'};
	die unless exists $data->{$model2}{'logL'};
	die unless exists $data->{$model1}{'logL'};
	die unless exists $data->{$model1}{'cnt_free'};
	my ($best, $diff, $chi) = chi_squared_test(	$data->{$model2}{'logL'},
												$data->{$model1}{'logL'},
												$data->{$model2}{'cnt_free'} -
												$data->{$model1}{'cnt_free'});

	if ($data->{$model2}{'prior_w'}[-1] < 1.0)
	{
		my $msg = ('8' x 60).
					"\nPositive test significance failed:\n\t".
					"Model M$model2 is best fit with an omega of < 1.0 ".
					"($data->{$model2}{'prior_w'}[-1])\n".
					('8' x 60). "\n";
		print STDERR $msg  if $verbose;
		print $output_file $msg;
		return 0;
	}

	if ($best != 1 || $chi == 0)
	{
		my $msg = ('8' x 60).
					"\nPositive test significance failed:\n\t".
					"Model M$model2 (logL = $data->{$model2}{'logL'}) ".
					"is not significantly better than ".
					"Model M$model1 (logL = $data->{$model1}{'logL'})\n".
					('8' x 60). "\n";
		print STDERR $msg if $verbose;
		print $output_file $msg;
		return 0;
	}
	my $msg = 	"Models M$model2 and M$model1:".
					"\tp < $chi".
					"\t(chi-squared diff = $diff)\n";
	print STDERR "\t", $msg if $verbose;
	print $output_file $msg;
	return 1;
}
#
#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888



#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#	test parameters for positive selection
#
sub test_parameters(\%$)
{
	my ($best_model_data, $PARAMETER_OUTPUT) = @_;

	my $success = 0;
	$success += test_model($PARAMETER_OUTPUT, %$best_model_data, 1,2);
	$success += test_model($PARAMETER_OUTPUT, %$best_model_data, 7,8);

	if ($success == 2)
	{
		print STDERR "\n    ", "8" x 30, "\n\n\tPositive selection!!!\n\n    ", "8" x 30, "\n\n" if $verbose;
	}
	else
	{
		exit(0);
	}
}
#
#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#	print parameters
#
sub print_parameters(\%\%$)
{
	my ($best_model_data, $model_data, $PARAMETER_OUTPUT) = @_;


	for my $model(1,2,7,8)
	{
		# save best logL file, then the rest
		my @files = ($best_model_data->{$model}{'file'});
		for my $file (keys %{$model_data->{$model}})
		{
			push(@files, $file) unless ($file eq $best_model_data->{$model}{'file'});
		}

		print $PARAMETER_OUTPUT "8" x 80, "\n>model M$model\n",
								"\tFile with best log likelihood printed first\n";
		print $PARAMETER_OUTPUT "8" x 80, "\n";
		for my $file (@files)
		{
			print $PARAMETER_OUTPUT ">model M$model from $file\n";
			print $PARAMETER_OUTPUT "\t Log likelihood  =\t", $model_data->{$model}{$file}{'logL'}, "\n";
			print $PARAMETER_OUTPUT "\t Free parameters =\t", $model_data->{$model}{$file}{'cnt_free'}, "\n";
			print $PARAMETER_OUTPUT "\t kappa           =\t", $model_data->{$model}{$file}{'kappa'}, "\n";
			print $PARAMETER_OUTPUT "\t tree length     =\t", $model_data->{$model}{$file}{'treelen'}, "\n";
			print $PARAMETER_OUTPUT "\t omegas          =\t", join (",", @{$model_data->{$model}{$file}{'prior_w'}}), "\n";
			print $PARAMETER_OUTPUT "\t probabilities   =\t", join (",", @{$model_data->{$model}{$file}{'prior_p'}}), "\n";
			if ($model_data->{$model}{$file}{'parameters'})
			{
				print $PARAMETER_OUTPUT "\t model parameters=\t", join ("; ", @{$model_data->{$model}{$file}{'parameters'}}), "\n";
			}
		}
	}
	print $PARAMETER_OUTPUT "8" x 80, "\n\n";
}
#
#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888



#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#

#	chi_squared_test
#
#		test for significance given log likelihoods and degrees of freedom

sub chi_squared_test($$$)
{

	my ($lnL1, $lnL2, $df) = @_;
	my $diff;
	my $best;
	my $chi;
	
	#chi-squared table
	my @chi_table = ();
		@{$chi_table[1]}	= ('3.84146', '6.63490' ,'10.828');
		@{$chi_table[2]} = ('5.99147', '9.21034', '13.816');
		@{$chi_table[3]} = ('7.81473', '11.3449', '16.266');
		@{$chi_table[4]} = ('9.48773', '13.2767', '18.467'); 	

	
	#calculate log_likelihood_difference X2
	if ($lnL1 > $lnL2)
	{
		$diff = sprintf("%.2f", 2*($lnL1 - $lnL2));
		$best = 1;
	}
	else
	{
		$diff = sprintf("%.2f", 2*($lnL2 - $lnL1));
		$best = 2;
	}

	#for a given degree of freedom, check backwards through the chi-squared table
	if ($diff >$chi_table[$df][2])
	{
		$chi = '0.001';
	}
	elsif ($diff >$chi_table[$df][1])
	{
		$chi = '0.01';
	}
	elsif ($diff >$chi_table[$df][0])
	{
		$chi = '0.05';
	}
	else
	{
		$chi = '0';
	}
	
	#return best model (1 or 2), 2 x log-likelihood difference, chi-squared probability
	return ($best, $diff, $chi);

}
#
#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888




#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#	Retrieve sequence length and gaps
#
#		Todo: Can add cDNA translatation here if desired
#
#
sub handle_gaps($\$\@)
{
	my ($FILE, $cnt_lines, $alignments) = @_;
	my @all_codons;
	my @all_codons_after_gaps;
	while (<$FILE>)
	{
		$$cnt_lines++;
		chomp;
		my $line = $_;
		next unless length $line;
		
		next if ($line =~ /seed used/);
		if ($line =~ /Before deleting alignment gaps/)
		{
			while (<$FILE>)
			{
				$$cnt_lines++;
				chomp;
				$line = $_;
				next unless length $line;
				unless($line =~ /^\s*\S+\s+((?:[A-Z\-]{3}\s+)*[A-Z\-]{3})/)
				{
					last;
				}
				my @codons = split /\s+/, $1;
				last unless (@codons);
				for my $i(0..@codons-1)
				{
					$all_codons[$i] .=$codons[$i];
				}
			}
			#print join ("\n", @all_codons), "\n";
		}
		elsif ($line =~ /^\s*\S+\s+((?:[A-Z\-]{3}\s+)*[A-Z\-]{3})/)
		{
			my @codons = split /\s+/, $1;
			@$alignments = ();
			return (scalar @codons, scalar @codons);
		}
		else
		{
			die "Error:\n\tUnexpected line at line $$cnt_lines:\n\t[$line]";
		}
		if (!@all_codons)
		{
			die "Error:\n\tCould not find codons at line $$cnt_lines:\n\t[$line]";
		}
		if ($line =~ /After deleting gaps/)
		{
			while (<$FILE>)
			{
				$$cnt_lines++;
				chomp;
				$line = $_;
				next unless length $line;
				last unless ($line =~ /^\s*\S+\s+((?:[A-Z\-]{3}\s+)*[A-Z\-]{3})/);
				my @codons = split /\s+/, $1;
				last unless (@codons);
				for my $i(0..@codons-1)
				{
					$all_codons_after_gaps[$i] .=$codons[$i];
				}
			}
			#print join ("\n", @all_codons_after_gaps), "\n";
			last;
		}
		else
		{
			die "Error:\n\tUnexpected line at line $$cnt_lines:\n\t[$line]";
		}
		if (!@all_codons_after_gaps)
		{
			die "Error:\n\tCould not find ungapped codons at line $$cnt_lines:\n\t[$line]";
		}
	}
	my $seq = '';
	my $pos = 0;
	my $pos_gaps = 0;
	my $len_gaps = scalar @all_codons_after_gaps;
	my $len = scalar @all_codons;
	for $pos_gaps(0..$len_gaps-1)
	{
		while ($pos < $len)
		{
			if ($all_codons[$pos] eq $all_codons_after_gaps[$pos_gaps])
			{
				$seq .= 'A';
				last;
			}
			else
			{
				$seq .= '-';
				$pos++;
			}
		}
		$pos++;
	}
	$seq .= '-' x ($len - $pos);
	#print "$seq\n$pos_gaps, $len_gaps, $pos, $len\n";
	@$alignments = get_alignment_marks($seq);
	return scalar @all_codons, scalar @all_codons_after_gaps;
}
#
#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888





#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#	Get best model based on log likelihoods
#
sub retrieve_best_model(\%\%)
{
	my ($best_model_data, $model_data) = @_;
	for my $model(sort keys %$model_data)
	{
		my $best_file;
		my $logL = 0;
		#print "model=$model\n";
		for my $file(sort keys %{$model_data->{$model}})
		{
			if ($model_data->{$model}{$file}{'logL'} < $logL)
			{
				#print "$model_data->{$model}{$file}{'logL'} < $logL for model $model $file\n";
				$logL = $model_data->{$model}{$file}{'logL'};
				$best_file = $file;
			}
		}
		#print "model $model, best file = $best_file\n";
		$best_model_data->{$model} = $model_data->{$model}{$best_file};
		$best_model_data->{$model}{'file'} = $best_file;
		#print "\nbest_model_data->{$model}{'logL'}= $best_model_data->{$model}{'logL'}\n";
	}
}
#
#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888








#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#	print positive residues as FASTA
#
sub print_positive_residues(\%$)
{
	my ($best_model_data, $PARAMETER_OUTPUT) = @_;

	open (POSITIVE_RESIDUES, ">$positive_residues") or
		die "Error:\n\tCould not open the file for saving positively selected residues in FASTA format [$positive_residues]\n$!\n";
	for my $model(2,8)
	{
		print POSITIVE_RESIDUES  ">model M$model from file $best_model_data->{$model}{'file'}\n";
		for (@{$best_model_data->{$model}{'probabilities'}})
		{
			if ($_)
			{
				print POSITIVE_RESIDUES int($_ * 10) % 10;
			}
			else
			{
				print POSITIVE_RESIDUES  "-";
			}
		}
		print POSITIVE_RESIDUES  "\n";
	}

	my $len = @{$best_model_data->{2}{'probabilities'}};
	my @consensus90 = ('-') x $len;
	my @consensus95 = ('-') x $len;
	my @consensus99 = ('-') x $len;
	my ($cnt_90, $cnt_95, $cnt_99) = (0, 0, 0);
	for my $pos (0.. $len -1)
	{
		#
		# use sorted list of probabilities
		#
		my @prob = sort { $a <=> $b}	($best_model_data->{2}{'probabilities'}[$pos],
										 $best_model_data->{8}{'probabilities'}[$pos]);

		# only need to test for probability of least significant

		# if you want a more sophisticated scheme can implement here
		# e.g. > 0.9 for one and 0.5 for the other two
		if ($prob[0] >= 0.9)
		{
			$consensus90[$pos] = '*';
			$cnt_90++;
			if ($prob[0] >= 0.95)
			{
				$consensus95[$pos] = '*';
				$cnt_95++;
				if ($prob[0] >= 0.99)
				{
					$consensus99[$pos] = '*';
					$cnt_99++;
				}
			}
		}
	}
	print POSITIVE_RESIDUES  ">Consensus (both models):  >= 0.9 [$cnt_90 sites]\n";
	print POSITIVE_RESIDUES join ("", @consensus90), "\n";
	print POSITIVE_RESIDUES  ">Consensus (both models):  >= 0.95 [$cnt_95 sites]\n";
	print POSITIVE_RESIDUES join ("", @consensus95), "\n";
	print POSITIVE_RESIDUES  ">Consensus (both models):  >= 0.99 [$cnt_99 sites]\n";
	print POSITIVE_RESIDUES join ("", @consensus99), "\n";

	print $PARAMETER_OUTPUT "\t$cnt_90 positively selected residues predicted by both models p > 0.9\n";
	print $PARAMETER_OUTPUT "\t$cnt_95 positively selected residues predicted by both models p > 0.95\n";
	print $PARAMETER_OUTPUT "\t$cnt_99 positively selected residues predicted by both models p > 0.99\n\n";
	if ($verbose)
	{
		print STDERR "\t$cnt_90 positively selected residues predicted by both models p > 0.9\n";
		print STDERR "\t$cnt_95 positively selected residues predicted by both models p > 0.95\n";
		print STDERR "\t$cnt_99 positively selected residues predicted by both models p > 0.99\n\n";
	}
}
#
#
#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888



open (PARAMETER_OUTPUT, ">$parameter_output") or
	die "Error:\n\t".
				"Could not open the file for saving the codeml parameters ".
				"for the different models [$parameter_output]\n$!\n";



print STDERR "\n\n" if $verbose;

my @filelist = @ARGV; #list of files for processing
	
#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#For each file, open and read contents into memory, and process
my %model_data;
my $seq_len = 0;
my %pos_prob;
my $len;
my $len_sans_gaps,
my @alignment;
foreach my $file (@filelist)
{
	open FILE, $file or die "Error:\n\tCould not open [$file].\n$!\n";
	
	my $model = 0;
	#my $omega_saved_for_model0 = 0;
	$seq_len = 0;
	my $cnt_lines = 0;
	
	#
	# find gap positions
	#
	my ($this_len, $this_len_sans_gaps) = handle_gaps(*FILE, $cnt_lines, @alignment);
	if (!defined $len)
	{
		($len, $len_sans_gaps) =($this_len, $this_len_sans_gaps);
		
		print PARAMETER_OUTPUT "Aligned sequence lengths = $len\n";
		print PARAMETER_OUTPUT "Aligned sequence lengths with gaps removed = $len_sans_gaps"
						if $len != $len_sans_gaps;
		print PARAMETER_OUTPUT "\n\n";
		if ($verbose)
		{
			print STDERR "\tAligned sequence lengths = $len\n";
			print STDERR "\tAligned sequence lengths with gaps removed = $len_sans_gaps\n"
							if $len != $len_sans_gaps;
			print STDERR "\n";
		}
	}
	else
	{
		if ($len != $this_len)
		{
			die "Error:\n\tThe sequence lengths in file [$file] are different. ".
				"Are you sure you used the same sequences throughout?\n";
		}
		if ($len_sans_gaps != $this_len_sans_gaps)
		{
			die "Error:\n\tThe ungapped sequence lengths in file [$file] are different. ".
				"Are you sure you used the same sequences and parameters throughout?\n";
		}
	}
	
	#
	# main processing loop
	#
	while (<FILE>)
	{
		$cnt_lines++;
		chomp;
		my $line = $_;
		next unless length $line;


		
		#store length
		if ($line =~ /ns\s*=\s*\d+\s*ls\s*=\s*(\d+)/)
		{
			if (!$seq_len)
			{
				$seq_len = pos_in_aligned_seq($1, @alignment);
				if ($seq_len != $len)
				{
					die "Error:\n\tThe lengths of sequences ($len) in file [$file] are different ".
						"from the sequence length used in the analysis ($seq_len). This is a bug!!\n";
				}
			}
			
		}
		
		
		#store model number
		elsif ($line =~ /^Model\s(\d)\:/)
		{
			$model = $1;
			#%model_data = ();
			%pos_prob = ();
			
		}
		
		
		
		#store number of free parameters and Log Likelihood
		elsif ($line =~ /^lnL\(ntime.+np\:\s+(\d+)\)\:\s+([\-0-9.eE]+)/)
		{
			$model_data{$model}{$file}{'cnt_free'} = $1;
			$model_data{$model}{$file}{'logL'} = $2;
		}
		
		
		#
		# store kappa
		#
		elsif ($line =~ /^kappa.+\=\s+(\d+\.\d+)/)
		{
			my $n = sprintf("%.2f", $1);
			$model_data{$model}{$file}{'kappa'} = $n;
		}
		
		#
		# store tree length
		#
		# tree length =  10.58821
		elsif ($line =~ /^tree length\s+=\s+([eE0-9.]+)/)
		{
			$model_data{$model}{$file}{'treelen'} = $1;
		}
		
		#
		# special handling for omega in model 0
		#
		elsif ($model == 0				&&
			   #!$omega_saved_for_model0 &&
				$line =~ /dN & dS for each branch/)
		{
			$line = <FILE>;$line = <FILE>;$line = <FILE>;$line = <FILE>;
			$cnt_lines +=4;
			$line =~ /\s*\d+\.\.\d+\s+[eE0-9.]+\s+[eE0-9.]+\s+[eE0-9.]+\s+([eE0-9.]+)/;

			push(@{$model_data{$model}{$file}{'prior_w'}}, $1);
			push(@{$model_data{$model}{$file}{'prior_p'}}, 1);
			#$omega_saved_for_model0++;
		}
		
		#
		# store parameters for beta model
		#
		elsif ($line =~ /^\s*p=\s+([eE0-9.]+)\s+q=\s+([eE0-9.]+)/)
		{
			my $p = sprintf("%.3f", $1);
			my $q = sprintf("%.3f", $2);
			push(@{$model_data{$model}{$file}{'parameters'}}, "p = $p", "q = $q");
		}
		
		
		#
		# store parameters for beta + w model
		#
		elsif ($line =~ /^\s*p0=\s+([eE0-9.]+)\s+p=\s+([eE0-9.]+)\sq=\s+([eE0-9.]+)/)
		{
			#parameters for beta+w model
			my $p0 = sprintf("%.3f", $1);
			my $p = sprintf("%.3f", $2);
			my $q = sprintf("%.3f", $3);
			push(@{$model_data{$model}{$file}{'parameters'}}, "p0 = $p0", "p = $p", "q = $q");
		}
		elsif ($line =~ /^\s*\(p1=\s+([eE0-9.]+)\)\s+w=\s+([eE0-9.]+)/)
		{
			my $p1 = sprintf("%.3f", $1);
			my $w = sprintf("%.3f", $2);
			push(@{$model_data{$model}{$file}{'parameters'}}, "p1 = $p1", "w = $w");
		}
		
		#
		# store total prior probabilities and omegas
		#
		elsif ($line =~ /^p:\s*(.+)/)
		{
			#my @p = split /\s+/, $1;	# fields can run into each other if > 99
			my @p = $1 =~ /(?:\d+\.\d{0,5})/g;
			push(@{$model_data{$model}{$file}{'prior_p'}}, @p);
		}
		elsif ($line =~ /^w:\s*(.+)/)
		{
			#my @w = split /\s+/, $1;	# fields can run into each other if > 99
			my @w = $1 =~ /(?:\d+\.\d{0,5})/g;
			push(@{$model_data{$model}{$file}{'prior_w'}}, @w);
		}
		
		#
		# store positively selected sites
		#
		#elsif ($line =~ /^\s+(\d+)\s\w\s(\d\.\d+)/)
		#{
		#	my $pos = pos_in_aligned_seq($1 - 1, @alignment);
		#	$pos_prob{$pos} =$2;
		#}

		elsif ($line =~ /^Bayes Empirical Bayes \(BEB\) analysis/)
		{
			$line = <FILE>;
			$cnt_lines++;
			chomp $line;
			next unless $line =~ /^Positively selected sites/;
			$cnt_lines += 3;
			for (0..2)
			{
				$line = <FILE>;
			}

			while (<FILE>)
			{
				$cnt_lines++;
				chomp;
				$line = $_;
				last unless ($line =~ /^\s+(\d+)\s+\w\s+(\d\.\d+)/);
				my $pos = pos_in_aligned_seq($1 - 1, @alignment);
				$pos_prob{$pos} =$2;
			}
		}
		
		#
		# End of model data: save probabilities
		#
		elsif ($line =~ /^Time\sused/)
		{
			die "Error:\n\tLog likelihood not parsed for model M$model in $file after line $cnt_lines\n\t "
					unless defined $model_data{$model}{$file}{'logL'};
			die "Error:\n\tFree parameters not parsed for model M$model in $file after line $cnt_lines\n\t "
					unless defined $model_data{$model}{$file}{'cnt_free'};
			die "Error:\n\tKappa not parsed for model M$model in $file after line $cnt_lines\n\t "
					unless defined $model_data{$model}{$file}{'kappa'};
			die "Error:\n\tTree length not parsed for model M$model in $file after line $cnt_lines\n\t "
					unless defined $model_data{$model}{$file}{'treelen'};
					
					

			#
			# Save sequence probabilities for positive sites
			#
			if (%pos_prob)
			{
				my @seq = (0.0) x $seq_len;
				for my $pos(keys %pos_prob)
				{
					die "$pos > $seq_len" if $pos > $seq_len;
					$seq[$pos] = $pos_prob{$pos};
				}
				
				push(@{$model_data{$model}{$file}{'probabilities'}}, @seq);
			}
		}
	}
}
#
#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888












#
#	Get best model based on log likelihoods
#
my %best_model_data;
retrieve_best_model(%best_model_data, %model_data);

#
#	test whether best model has +ve selection
#
test_parameters(%best_model_data, *PARAMETER_OUTPUT);
print PARAMETER_OUTPUT "\n\n";


#
#	print parsed parameters to file
#
print_parameters(%best_model_data, %model_data, *PARAMETER_OUTPUT);


#
#	print positive residues as FASTA
#
print_positive_residues(%best_model_data, *PARAMETER_OUTPUT);


